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Fluctuations in systems away from thermal equilibrium have features that have no analog in 
equilibrium systems. One of such features concerns large rare excursions far from the stable state 
in the space of dynamical variables. For equilibrium systems, the most probable fluctuational 
trajectory to a given state is related to the fluctuation-free trajectory back to the stable state by 
time reversal. This is no longer true for nonequihbrium systems, where the pattern of the most 
probable trajectories generally displays singularities. Here we study how the singularities emerge as 
the system is driven away from equilibrium, and whether a driving strength threshold is required for 
their onset. Using a resonantly modulated oscillator as a model, we identify two distinct scenarios, 
depending on the speed of the optimal path in thermal equilibrium. If the position away from the 
stable state along the optimal path grows exponentially in time, the singularities emerge without 
a threshold. We find the scaling of the location of the singularities as a function of the control 
parameter. If the growth away from the stable state is faster than exponential, characterized by the 
ability to reach infinity in finite time, there is a threshold for the onset of singularities, which we 
study for the model. 



I. INTRODUCTION 

The phenomena of large rare fluctuations have been 
attracting increasing attention recently in the context of 
small systems, such as nanomechanical resonators [11 12], 
nanomagnets [3", 4 , Josephson junctions [5^, and trapped 
ions [6 . Another arena of rare events is biology. De- 
spite being inherently stochastic [7], biological processes 
are notoriously robust [8','9^, yet prone to mutations, and 
hence to disease, genetic variability, evolution, as well 
as making functional use of the noise [10] . Large fluctua- 
tions may also manifest themselves in catastrophic events 
in various phenomena, such as in climate [11 , or popu- 
lation dynamics [12 . 

The problem of fluctuations has two closely related 
aspects: the form of the probability distribution and 
properties of paths in the dynamical configuration space 
followed during a fluctuation. Both aspects are well- 
understood for systems in thermal equilibrium. Here the 
probability density distribution is given by the Boltz- 
mann law. When the fluctuations caused by thermal 
noise are weak on average, this distribution is narrowly 
peaked at the stable state of the system in phase space 
and usually rapidly falls down away from it. The paths 
followed in fluctuations were first studied by Onsager and 
Machlup [I3l [14], who showed that in a large fluctua- 
tion to a remote part of the phase space the system is 
most likely to follow a certain optimal path. For sys- 
tems in thermal equilibrium this path is simply related 
to the fluctuation- free path (may also be called relax- 
ational or mean-field path) back to the stable state via 
a time-reversal relationship. This relationship follows 
from the principle of maximal work for equilibrium sys- 
tems |15j. It implies that if the pattern of relaxational 
paths is unique, i.e. relaxation from a given point to 
the attractor is accomplished via a unique path (this is 
a generic situation), then the pattern of optimal fluc- 



tuational paths, without conditioning on the fluctuation 
time, will also be unique - a given point in phase space 
is reached from an attracting state by a unique fluctua- 
tional path. In contrast, in nonequilibrium systems there 
is no simple relation between the most probable fluctua- 
tional and relaxational paths. Moreover, generically the 
pattern of optimal fluctuational paths may allow for in- 
tersection of neighboring paths, and the loci of points 
in phase space along which this happens are the caustic 
singulatiries [16l |T7]. Such singularities are a counter- 
part of caustics in optics [18] and the WKB trajectories 
in Quantum Mechanics Caustics are accompanied 

by switching lines which separate regions of phase space 
such that the neighboring points on the opposite sides of 
switching lines will be reached by non-neighboring fluc- 
tuational paths. The gradient of the probability distribu- 
tion has discontinuities along switching lines. Switching 
lines are the experimentally-observable signatures of sin- 
gularities. The existence of singularities implies a stark 
qualitative difference from equilibrium - both in the pat- 
tern of optimal paths and in the form of the probability 
distribution. 

Such picture of the structure of singularities is under- 
stood [17^. The separate question of the mechanism by 
which they are formed has been studied for systems with 
a saddle point by Dykman, Millonas, and Smelyanskiy 
[TT] (see also [20]) and by Maier and Stein (for example 
in [2r , among others) , and for systems with an unstable 
limit cycle in [22, . Here we focus specifically on the ques- 
tion of threshold in a monostable system as it is driven 
away from equilibrium. We frame the study of singular- 
ities of optimal paths in the following context. Consider 
a generic system with a finite-dimensional phase space. 
Driving of a system from equilibrium will be character- 
ized by a certain control parameter that serves as a quan- 
titative measure of the deviation from equilibrium. Sin- 
gularities will generically appear at some critical value of 
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this control parameter. With this picture in mind, we 
formulate our main question as follows: what properties 
of the system determine whether this critical strength 
of the driving is finite or infinitesimal? The effect of 
the perturbation to first-order is computed along the un- 
perturbed path, and it is given by a certain functional 
of the unperturbed path. Therefore, the answer to our 
main question depends on the form of the perturbation 
and on the properties of unperturbed paths. Focusing on 
the case of non-singular perturbations, we point out the 
qualitative difference between those unperturbed paths 
that reach infinity in an infinite or finite time, respec- 
tively giving rise to a zero versus a finite critical driving 
threshold at which singularities appear. 

In this paper we study how singularities in the pat- 
tern of optimal paths emerge in a particular system as 
it is driven away from thermal equilibrium - a nonlinear 
oscillator driven by periodic field and interacting with a 
thermal bath. The amplitude of the periodic driving field 
is a control parameter that serves as a quantitative mea- 
sure of the deviation from equilibrium. The situation at 
zero field is straightforward and agrees with predictions of 
equilibrium statistical mechanics: the probability distri- 
bution in the phase space of the oscillator is predictably 
Boltzman and a most-likely fiuctuation to any specified 
point of phase space is realized via a unique path. More- 
over, these paths can be obtained from the noise- free 
paths by reversing the time of the non-Hamiltonian part 
of the dynamics. In contrast, the situation at finite driv- 
ing is far less clear. We focus on the question of the 
driving threshold above which the singularities develop 
- is it finite or infinitesimal, and what properties of the 
system control it? 

In our oscillator case-study, the presence of finite or in- 
finitesimal driving threshold turns out to depend on the 
strength of the nonlinear damping. When the nonlinear 
damping is present, the value of the threshold driving 
field required for the appearance of caustics is nonzero. 
This threshold field scales linearly with the strength of 
the nonlinear damping and becomes zero when only the 
linear damping is present. In this regime of purely lin- 
ear damping, caustics appear at an infinitesimal value 
of the driving field. This is a very interesting scenario: 
the system becomes qualitatively different from an equi- 
librium one even if it is driven away from equilibrium 
by an infinitesimal amount - a type of an equilibrium 
fragility. Secondly, we found that as the driving field 
passes above the threshold, the caustics arrive from in- 
finity and move in closer to the attracting stable state as 
the driving field increases. We used perturbation theory 
to understand the mechanism of caustic formation. We 
learned that caustics form as a result of the skewness (ra- 
dial asymetry) of the pattern of fluctuational paths. This 
skewness decreases with decreasing driving field, and this 
is why caustics are found further away from the stable 
state as driving goes to zero. The analytical treatment 
also explains the appearance of a finite threshold when 
nonlinear damping is finite. The presence of nonlinear 



damping causes faster than exponential dynamics and 
therefore allows paths to reach infinity in finite time. 
The interplay between the skewness, which tends to cause 
caustics, and this accelerated motion in the presence of 
nonlinear damping is responsible for the threshold. We 
emphasize that the ability to reach infinity in finite or in- 
finite time serves merely as a way to separate the speed of 
optimal paths into two distinct categories. The physics 
does not depend on the point at infinity, but it does de- 
pend on the speed of the optimal paths in the relevant 
parts of phase space. In fact, the model needs to be 
modified far away from the stable state if it is to remain 
realistic. However, as the driving field grows, the caustics 
are situated sufficiently close to the stable state where the 
present model is realistic. 

This paper is organized as follows. Section|ll|is devoted 
to the exposition of the phenomenology of singular fea- 



tures. In Section [ni| we outline our analytical approach 
and compare it with numerical predictions. We further 
discuss the relevance of this work in Section HVl 



II. 



PATTERN OF FLUCTUATIONAL PATHS IN 
THE PHASE SPACE OF A NONLINEAR 
OSCILLATOR 

A. Linear Damping 



We consider a classical nonlinear oscillator coupled to 
a heat bath 
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aq 



F cos (Qt) ^ f{t) 



(1) 



where the effect of the bath has been reduced to the 
damping and noise terms which satisfy a fluctuation- 
dissipation relation 



{f{t)f{t'))=2jkT5{t-t'). 



(2) 



This description can be obtained, for example, by treat- 
ing the bath as a collection of harmonic oscillators, each 
coupled linearly to the oscillator described by the vari- 
able q [23^. In the regime when the resonance is sharp: 
7/cJo <C 1 and the oscillator is driven close to resonance: 
|1] — cjol ~ 7, its dynamics is given approximately by fast 
simple harmonic motion modulated by a slow function of 
time, which is described by the following amplitude equa- 
tions in "slow time" T (see Appendix A) : 
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(/,(T)/,(T'))=D5,,,(5(T-T'), 



(8) 



which can ah be obtained from the microscopic Hamilto- 
nian by integrating out the bath degrees of freedom [23 
or from the corresponding Langevin equation, Eq. (IlJ 
[24-26 . The functional probabihty density of a given 
reahzation of a fluctuation can be expressed as 



p{path) ~ exp 



s [fQ{t),fp{t), \Qmpit),Q{t), p{t)] \ 



D 



(9) 

For Gaussian noise with the correlator {fi{t)fj{t')) = 
DdijS{t — t') (white noise), the functional S is given to 
lowest order in D by [27H29] : 



3 = {Q.P} 
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iXjit) {xj - Kj - fj{t)) 



dt 



(10) 

where x stands for either Q{t) or P{t) and Ks are also 
functions of these dynamical variables. Performing the 
functional integral in As and /s leaving the following 
Wentzel-Freidlin action 130 . 



S 



dt. (11) 



Because our goal in this work is to understand the prop- 
erties of the exponential part of distributions in the limit 
of weak noise strength D we follow the WKB approach - 
evaluate the functional along the paths which minimize 
5, or optimal paths. The problem then boils down do to 
finding optimal Q{t) and P{t) which minimize the asso- 
ciated action S. Interpreting the integrand in Eq. (11) as 



the Lagrangian, the associated Wentzel-Freidlin auxiliary 
Hamiltonian is 



(12) 



and the optimal paths are solutions to the correspond- 
ing Hamiltonian motion generated by this auxiliary 1-L. 
We will now study the properties of these trajectories 
numerically, and then analytically, with the help of the 
perturbation theory in Section III A 



1 . Numerical procedure 

We note that p = ^{x—K) = noise, so in the absence of 
noise, the only possible dynamics lie on the p = plane 
- the "mean- field" relaxational dynamics. Any fluctu- 
ational dynamics will be lifted out of this plane. The 
fixed points (FP) of this dynamics are (Qfp, ^fp, 0, 0) 
where {Qfp^Pfp) are the FP of the original 2D dynam- 
ics. It can also be easily shown that eigenvalues of each of 
these auxiliary FP are (Ai, A2, — Ai, — A2) where Ai and 
A2 are eigenvalues of the FP of the original dynamics. 
The p = plane forms the stable manifold of each of the 
fixed points. A trajectory that escapes from an attract- 
ing fixed point will lie on a certain manifold in the auxil- 
iary (Q, P^pq^pp) space. In the vicinity of the attracting 




FIG. 1. An envelope of intersecting neighboring straight lines 
is a simple example of caustic, as seen at the top of the figure. 



fixed point, this manifold is tangent to the plane spanned 
by two unstable eigenvectors of this FP. Thus, the fluctu- 
ational trajectories lie on the unstable "Lagrangian Man- 
ifold" j3TH34] of the attracting fixed point. We calculate 
the set of trajectories that escape the attractor by evolv- 
ing a set of initial conditions placed on a circle on this 
unstable manifold very close to the attracting FP. 

All calculations of trajectories are terminated either 
when a stopping radius is reached or when a trajectory 
reaches a caustic. In section III A we show explicitly 



that radius of a trajectory at = grows as a simple 
exponential, which means that to provide an accurate 
solution valid at larger radii, the time steps must de- 
crease dramatically. The choice of the stopping radius 
was dictated mostly by the position of interesting fea- 
tures, such as caustics, if they are present; otherwise, 
such as in Fig. [2]- (a), a convenient stoping radius was 
chosen. For a given value of stopping radius, the size of 
time step was decreased until the results did not change 
significantly. Fourth-order Runge-Kutta numerical inte- 
grator was used. 

A caustic is a the set of points where neighboring 
trajectories intersect - Fig. [l] In our case we are con- 
cerned with the intersection of projections of trajectories 
of the auxiliary Hamiltonian dynamics. Intersection of 
trajectories are of course an artifact of the projection - 
the evolution in the 4-dimensional auxiliary phase space 
with coordinates {Q-,P-,Pq-,Pp) is smooth and unique. As 
we mentioned, the trajectories in question evolve on the 
2-dimensional Lagrangian manifold, which is a surface 
p(r), i.e., (pq(Q, P),Pp(Q, P)), which emanates from the 
fixed point. In the vicinity of the fixed point, this man- 
ifold is a plane, but further away it may have folds. As 
two neighboring trajectories pass over this fold, their pro- 
jections will intersect; excellent cartoons of this can be 
found in [TZ and also in the beginning of ^35^. The fold 
is characterized by such location where the tangent plane 
to the surface is perpendicular to the (Q, P) plane. This 
happens when 



J d{pQ,pp) 



(13) 



SO caustics are mapped out from the knowledge of p(r). 

The function p(r) is not known over all r - only p along 
a given trajectory can be calculated. However, as the 
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equations of motion for pg, pp, Q and P are integrated 
in time, it is also possible to simultaneously integrate 
the equations of motion for second derivatives such as 
QQQP = ^§§-1 etc. These the time-dependent Riccati 
equations, [32j. 



dpkdpi 



dvjdpk 



dvidpk 



dvidrj 



(14) 



With the help of these and the other four dynamical 
equations of motion we can track the evolution of J: all 
together, there are 7 coupled ordinary differential equa- 
tions. The initial conditions are obtained from the knowl- 
edge of the surface p(r) in the vicinity of the attractor, 
defined by unstable eigenvectors. 



2. Optimal fluctuational trajectories - numerical results 

In this subsection we will study the pattern of fiuc- 
tuational trajectories at ever increasing driving strength 
T to elucidate the physics. When = 0, the pattern 
is rotationally symmetric. Fig. [2]- (a), a fact that will be 
proven below. As becomes non-zero, the rotational 
symmetry of the pattern of trajectories becomes broken. 
The density of trajectories also no longer obeys radial 
symmetry. If an observation is made at a given radius, 
one would notice an ever increasing density of trajecto- 
ries in certain regions until they actually "jam" into each 
other, and subsequently intersect. This scenario is de- 
picted in Fig. |2]-(b), along with the numerically mapped 
location of the caustic, obtained by the method described 
Section |II A 1| The caustic displays a cusp structure, 
which is expected on the generic basis (see [17], the Ap- 
pendix 12 of |34| and [36l[^). As is increased further, 
the cusp continues to pull in closer to the origin and also 
becomes more open. We summarize the relationship be- 
tween the radius of the cusp and in Fig. |3] for two 
values of 7^. In this log- log plot we also show that each 
of the numerical sets of data is well-approximated by a 
straight line with slope —1/2, a power law that we predict 
analytically in section [ill A| As T continues to increase, 
the caustics evolve to become more complex. We will 
not focus on the details of this in current work - our goal 
here is to focus on those features which are expected to 
be generic, such as the scaling properties of singularities 
as ^ and to answer the question whether there is a 
finite threshold in T at which the caustics appear. 

Varying r] at fixed T also reveals an asymptotic scaling. 
Fig. [4j For the value of shown, the relationship of R 
versus r] becomes asymptotic to a power law for larger r^. 



= 1 + Ar]P, 



(15) 



where the exponent p ^ 1.26 for this value of J^. More- 
over, since ^ = where 6uj is the deviation of the driv- 
ing frequency from the resonant frequency, we also learn 
that R (Suj) is asymptotic to a power law near resonance. 
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FIG. 2. (a): The pattern of 10 escape trajectories at = 0, 
rj = 0.4. The slow-down radius is 1. (b): The pattern of 10 
escape trajectories at = 0.008, r] = 0.4. Increased density 
of trajectories is clearly seen. The trajectories intersect along 
a caustic, shown in red, which originates at the cusp. Due to 
the spiral nature of the trajectories, their intersections can not 
be easily seen by eye - the caustic has been mapped out with 

(b)-inset: 



the help of the technique outlined in Secion II A 1 



As is increased the caustic system evolves and the cusp 
becomes more open. Here the small portion of the caustic 
around the cusp point is shown for = 0.05, ?7 = 0.4. At 
this value of the caustic system becomes very complicated. 
The full details of it are not shown, and are not expected to 
be generic; only the cusp part of the caustic is depicted. 
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FIG. 3. Log-log plot representing the radius of the cusp, 
measured from the origin versus the strength of the driving 
for ?7 = 0.4 circles and 77 = 0.1 diamonds. The dashed 
lines represent power laws with exponent —1/2, a theoretical 
prediction outlined in the next section. 




FIG. 4. Radius of the cusp versus 77 for a fixed 0.0012. 
The asymptotic exponent p in = 1 + Ar/^ is ~ 1.26. 



B. Nonlinear Damping 

We can also consider the nonlinear interaction with the 
bath. The Langevin description of this process is 

q^-fq^^q^q^ulq^aq^ = F cos (m) + + (16) 

with 

{mf{t'))=2jkTSit-t'). (17) 
The corresponding slow-time amplitude equations are 



dQ 
dT 
dP 
dT 



Kq[Q,P) + fQ(T) + QfQ{T)-Pfp{T) (18) 
Kp{Q,P) + fp{T) + P~fQ{T) + QfpiT) (19) 



Kp = 



dg 
dQ 

{MT)fjiT'))=DS,jSiT-T'), 



(20) 
(21) 
(22) 



/i is defined in the Appendix, and D/D = r]/fi. The ac- 
tion S now contains new nonlinear damping and noise 
terms associated with nonlinear interaction of the oscil- 
lator with the bath 

1 /"OO -| POO 

S= j {fUt) + fl(t)) dt+j} / (flit) + flit)) dt 

- r i^Q(t) (Q-Kq- fgit) - QfQit) - Pfp(t)) dt 

J —00 

- r i\p{t) (P-Kp- Qfp{t) + PfQit)) dt + 0{D) 

J — 00 

(23) 

The corresponding Wentzel-Preidlin Hamiltonian is 



PqKq +ppKp. 

(24) 



1. Optimal fluctuational trajectories - numerical results 

We now present our numerical findings pertaining to 
the case of non-zero ja. The key observation is the ex- 
tinguishing of caustics at large-enough values of ja. This 
is depicted in Fig. |5) where we plot the inverse of the 
radius of caustic versus fi for two different values of T. 
The J^-dependence at zero // was depicted in Fig. [3] The 
data is very well approximated by a fit of the form 



Rcaustic ^ {l^thr — 



-1/2 



(25) 



strongly suggesting the existence of a critical value of 
nonlinear damping jithr {^) , above which the caustics are 
pushed away to infinity. The existence of the such criti- 
cal fi can also be inferred from the trend of the data of 
Rcausitc versus /i, as seen in Fig. [H] The quantity (ithr 
versus is depicted in Fig. |6] for low values of /i and J^. 
For a given value of there is a critical value of /i below 
which caustics exist, iithri^)- This jithr appears to scale 
linearly with at low T 



lim /ithr{J^) 



(26) 



Equivalent ly, for a given value of /i, a finite value of is 
required for the formation of singularities, J^thr{l~^)' This, 
Tthr appears to scale linearly with ji at low ji: 



lim Tthri/^) ^ 
This fact combined with Eq. ( [25| ) implies 

^caustic ("^ thr^ 



(27) 
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of the field is manifested in the fact that the caustic itself 
is pushed into far regions of the phase space. However, 
quantities which obey power laws, such as the charac- 
teristic location of the caustic versus the strength of the 
driving field are nevertheless amenable to perturbative 
estimates. In doing this, we extract the characteristic 
location of the caustic indirectly, by identifying those re- 
gions of the phase space where the zeroth-order term in 
the action becomes comparable to the first-order correc- 
tion due to the driving. Using this technique we were 
able to reproduce several scaling laws with respect to the 
driving field and with respect to the strength of the non- 
linear damping. 



FIG. 5. Inverse of the caustic radius R versus the strength 
of the nonlinear damping /i for two different values of the 
driving force J^. Plotting 1/R helps to see the existence of 
the limiting value of /a above which caustics disappear. Solid 
curve - fit to the form {/j^thr — /x)"*^^^. 




FIG. 6. Numerically-obtained /Jthr verus (dots). At a 
given the caustics appear for /i < /j^thr- Alternatively, at a 
given /i, the caustics appear for > Tthr- The error bars are 
obtained from the uncertainty of exact value of the threshold, 
as seen in Fig. [5] The dotted line has a slope 6.4. 



III. ANALYTICAL SCALING ARGUMENTS 

The goal of the analytical treatment is to (i) explain 
various results of the previous section, but more im- 
portantly (ii) develop an intuitive understanding of the 
mechanism of the development of singularities. This 
should serve our ultimate goal - to understand what 
properties of the system control the singularity thresh- 
old value of the tuning parameter that takes the system 
away from equilibrium. The study of optimal paths will 
be based on the perturbation around the optimal path 
at zero J^, which we know from general arguments of 
statistical physics and can also calculate exactly. The 
perturbation theory can not be used directly to identify 
the location of the caustic because the effect of the per- 
turbation is always large on the caustic [16 ; the weakness 



A. Linear Damping 

In the limit of small noise the optimal path that 
minimizes S dominates the functional integral, and a 
quadratic-order WKB approximation may be used (with 
some exceptions, such as near caustics [ISj)- Applying 
the functional derivative to we generate the following 
set of Euler-Lagrange equations for the optimal trajec- 
tory 
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(29) 
(30) 
(31) 
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(33) 
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AO) 
JQ 
AO) 
Jp 



dP dP 
Using this it is easy to prove that at zero driving the 
optimal noise paths satisfy 

from which it would also follow that 

dQ 

This is in accord with the general "Principle of Max- 
imal Work" [15 for equilibrium systems - the optimal 
paths obey equations where the non-Hamiltonian part of 
the noise-free dynamics has been time-reversed; switch- 
ing the sign of rj accomplishes the same task as perform- 
ing the following transformation: T — T, Q ^ Q and 
P — P, but his is just a time reversal relationship, as 
can be verified from the original Duffing equation ([T]). 



-vQ 

rjP. 



(35) 
(36) 

(37) 
(38) 



for the action reduces to 

1 

4 



In view of Eqs. (10) and (29)-(30), the full expression 



S= - 



ifQit) + fkt)) dt, 



(39) 
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although to evaluate this integral we need a full expres- 
sion for /q and /p, which we do not know apriori. We 



note however, that at zeroth-order we can use Eq. (35) 
and Eq. (36) to obtain 



Sio) 



V 



(r(°)(i))^ 



dt 



(40) 



The solution Eqs. ^ -g for P^°\t) and Q^°\t) is 
most conveniently expressed in polar variables: 



(P^°\t) = (Po + t - R 



2 i^' "' - 1) 
27? 



(41) 
(42) 



Combining this result for A^\t) with Eq. (40) we get 



5(°)(r,^) = ^ 



(43) 



where we switched notation from (i^, ^o) to (r, 6>) - such 
switch from treating a point in phase space as an initial 
condition to treating it as an independent variable will 
prove useful again below. We mention in passing that 
one can also use the Wentzel-Freidlin action, Eq. (11) as 



the starting form for the analysis. For example, the time- 
reversal symmetry of the fluctuational and relaxational 
paths can also be proven using this form as the starting 
point ^24j . 



1. Analysis of the first- order perturbation 

Since first order variation of L^^^ around the unper- 
turbed path is zero by definition [15 , the only possible 
first order in correction to the action is 



5^=/ L(i)(Q(°)(t),p(°)(t),...) 



dt. 



(44) 



The first-order correction to the action is the integral 
of the first order correction to the Lagrangian along the 
unperturbed path. In view of this and of Eqs. ([5|-(|6| and 
Eq. (10) we have 



5 = 5(°)+^5« + ... 
SW = f° ^pio){t)dt 

J — oo 



(45) 
(46) 



where we also made use of Eq. ( 32 ) and Eq. ( 36 ). We now 



analyze this contribution. With the help of Eqs. (41 )-(42 ) 
we have 



rjre'''' sm I t/ + t (e^"^^ - l) ] dt 

-oo V ) 

(47) 

Notice that the dependence on the field point (r, Q) has 
been put into the parameters of the integrand, allowing 
the limits of integration to be simple. We can change 



variables to x = ine^^^ and re-express the form of S^^^ in 
a more compact form: 
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{^f^y Jo 
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(48) 
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where ^ = ^ + ^ and = When r ^ y/rj the upper 
limit can be replaced by oo (notice that z is independent 
of r, so as r is increased but r] is held fixed, only the 
upper limit changes, but the structure of the integrand 
does not, justifying this replacement). The integral is 
taken in the x-plane along the imaginary axis, but the 
integrand is an analytic function in the right half-plane, 
so we may rotate the contour to go along the real x-axis. 
Then the integral coincides with the definition of the F- 
function in the variable z. After substituting hz and a bit 
of algebra we have. 



^(1) ^ V^^\'+l 



e^.r(^ + ll+c.c. 



(49) 

We note that limit ^ oo is actually equivalent to start- 
ing out with the integral that does not contain the In- 
term. The S'^^^ just obtained is a non-growing function, 
but with ever increasing oscillation frequency in the ra- 
dial variable. Next, we take a Jacobean, Eq. ([Tst. Using; 

Jacobean a compact form 



the series expansion for 5, and Eq7T43|), the 



After some algebra we see that at large radii. 



1/2 



X g4r7 



2r] 



(50) 



(51) 



Recall that exact J must diverge on caustics. In the 
language of perturbation theory, however, J is unlikely 
to diverge (clearly the perturbative solution to the action 
S'^^^ + S'^^^ is smooth and can not contain singularities; 
in fact, the Gaussian WKB theory is not meant to work 
on caustics [16]). Therefore, the signature of caustics 
in the language of first-order perturbation theory is not 
the divergence of J, but the point at which such theory 
begins to break down. In other words, the signature of 
caustics is the condition ~ l^^"*^^]- Thus, Eq. (51) 

suggests that 



5/4 



jri/2 



X e 8^ 



2n 



-1/2 



(52) 



In view of Stirling's approximation the last factor ap- 
proaches (27r)~^/^ in the limit of 77 ^ and 7r~^/^ in the 
limit ?7 ^ 00. Numerical computations show that this 
factor remains very close to the former constant value 
until 7? ~ 1, after which there is a crossover into the lat- 
ter value. Therefore the square root factor in Eq. (52) 



can be disregarded, and we recover the numerical result 
for scaling with respect to the driving field and friction. 



R, 



caustic 



n 



,5/4 



jri/2- 



(53) 
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B. Nonlinear Damping 



The analysis of II A can be repeated for the case of 
nonlinear damping. As in the linear case, we can take 
variational derivatives (there are two more now, with re- 
spect to the new noises), and find optimal /(t)s, /(t)s, 
A(t)s and {Q{t),P{t)). First, we find that at J* = 



= 4/iQP. 



(54) 
(55) 



in addition to Eqs. (35)- (36), which once again imply 
that the optimal paths at zero driving are time-reversed 
relaxations. Along the optimal path the full action is 
given by 



S 



) ^ M J — oo ^ ^ 

(56) 

We can also find zeroth-order trajectories: 



R 



Combining all these results we get 
same as in the linear damping case 



V 



(1 



(57) 



(58) 



(59) 



1. Analysis of the first- order perturbation 

We move on to calculate S'^^^ for the case with nonlin- 
ear damping. The statement in Eq. ( |44l) r emains true. 
Then, in view of Eqs. (18)-(19) and Eq. (|23p we have, as 
in the linear case 



^(1) = f° ^p(0)^t)dt, 

J — oo 



where we also made use of Eq. (32) and Eq. (36). With 



the help of Eqs. (57)-(58) and changing variables to ^ = 



5(1) 



re 
~2i 



1 ' ixr^ 



f/" 1 + ^ (1 - e') 

V ^ 



+ c.c. 
(60) 

This is the general form. We simplify notation by in- 

2 

troducing a = and A = 1 + a. There are now two 
important regimes to consider: regime where the nonlin- 
ear damping dominates the linear damping, a 1, so 
^ ^ 1 — ^ ^ 1 and the regime of perturbative non- 
linear damping, where it is weak compared to the linear 
damping, a<Cl, so^~a<Cl. 



In the first regime, the second term in the brackets in 
Eq. (60) dominates. This allows all the r-dependence to 



be taken out of the integral, giving 



5(1) 



r 



(61) 



We apply Eq. (|50|) to take the Jacobean, and find that 

-2 



r 



(62) 



This explains the existence of a threshold in J^. A suffi- 
ciently large is necessary at any large r to bring J^J^^^ 
to be comparable to J*^^^ {= rf). In contrast, the case of 
linear damping taught us that J^^^ is an always increas- 
ing function (~ r^), so for any finite value of there 
exists a radius where FJ^^^ dominates J^^\ 

We now consider the second mentioned regime. For 
this task, we rewrite Eq. (60) as follows 



re 



2i 



c.c. 



(63) 

Analogously to the case with linear damping, now define 
a new variable y = where = Also notice that 
in the limit a ^ 0, we have a/ A a. Then, 

JO 



4i Jo 



(1 - 2/i?/)2M 2 dy^c.c. 

(64) 



As in the case of linear damping we assume: r ^ y^, 
i.e. hi ^ 1. The integral in Eq. (64) is then a function 



of only /i and r^, so it will not affect our quest for scaling 
with r. Thus, 



5(1) 



where 



■In 1+^ 



1 



2r] J 



lim /(/i,7^) 



2ri ^ 

g4r7g 4 r 



2r} 



In 



'"^/(/^,^)+C.C. 

(65) 



c.c. 



(66) 



This S^^^ agrees with Eq. (49) as ji 
lirrie^o (1 — ex)^^"^ 



(the identity 



may be helpful). Expanding 
the logarithm in powers of /i and doing more algebra 
we find that a /i-contribution to the Jacobean with the 
largest power in r is 



J'/i^ X e 



/o(^) 



(67) 



This is the leading-order addition to the linear damping 

2 

Jacobean when r is large (i.e. r ^ V^)' while ^ is 
small. In other words, as we go to large r and decrease fi 
to be in the regime of perturbative effect of the nonlinear 
damping, the Jacobean will take the following asymptotic 
form [c.f. Eq. (IsTl)] 



.1/2 



,3/2 



"2^ 



(68) 
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Comparing with the magnitude of the J^-term now 
gives (omitting factors of r]) 



1 



11 



caustic 



(69) 



Nonhnear damping pushes caustic to larger radius. This 
describes the small- /i part of Fig. [5] In fact, combin- 
ing Eqs. (25) - (27) and expanding it in series in /i we 
reproduce this result. However, the perturbation the- 
ory does not give prediction for any value of /i, such as 
empirically-obtained Eq. (25). 



IV. SUMMARY AND DISCUSSION 

We considered a simple system interacting with a heat 
bath. The amplitude of the periodic driving is the 
control parameter that measures the deviation of the sys- 
tem from thermal equilibrium. We may summarise our 
findings as follows. At = the pattern of optimal fluc- 
tuational paths is smooth, and is related to the pattern 
of relaxational paths via the change of the signs of lin- 
ear and nonlinear damping. A unique fluctuational path 
leads to a given point (Q, P) in the phase space in the 
stationary regime. In the case of linear damping, the op- 
timal paths at = have exponentially-growing radii, 
and hence reach infinity in infinite time. With finite non- 
linear damping /i, the exponential growth is replaced by 
a faster growth beyond a certain radius, such that the 
optimal path reaches infinity in a finite time. This quali- 
tative difference, already present at = 0, has the effect 
upon the presence of finite or infinitesimal Tthr necessary 
for the appearance of singularities. 

In the case of /i = 0, we observe the appearance of sin- 
gularities for an infinitesimal value of The scenario by 
which this happens is such that the singularities appear 
at infinity and move in to smaller radii as is increased: 
Rcaustic ^ T~^l^. Interestingly, the singularities hap- 
pen not where the rotation switches from clockwise to 
conter-clockwise motion, see Fig. |2j after all, the region 
where the angular part of the dynamics slows down is ex- 
pected to be most sensitive to perturbations. We learned 
that the mechanism for the formation of singularities is 
instead related to skewness, or uneven distortion of the 
family of optimal paths. The smaller the T ^ the smaller 
is this distortion, and the larger is the distance required 
for the distortion to cause neighboring paths to intersect. 
To verify that the slowing-down region is non-essential for 
the appearance of caustics, we considered a system where 
g is given by 



1 



QT 



(70) 



(compare to Eq. ([t])). The optimal paths of this system 
rotate only in one direction and do not display a turn- 
around effect as do the paths we have studied so far. 
However, we found the caustics in this system as well. 



In the case when the nonlinear damping /i is non-zero, 
we found that 



Rcaustic 
^ thr 



thr) 



-1/2 



(71) 
(72) 



Thus, for a fixed value of /i, varying T has the effect 
of bringing the caustics in from infinity. We can under- 
stand the origin of the linear relationship Tthr ^ M from 
the following argument. There are two effects. On the 
one hand, a non-zero driving T introduces skewness of 
the pattern of optimal paths. The characteristic radius 
at which the skewness is sufficient to create a singularity 
scales like T~^l^ at zero /i. On the other hand, a finite 
value of \i causes the optimal paths to accelerate; their 
dynamics switches from an exponentially growing to an 
even faster one, causing them to reach infinity in finite 
time. This change in the dynamics of optimal paths hap- 

1/2 

pens at a characteristic radius ~ (^/m) • When the 
motion is faster than exponential at the radius where 
the skewness would otherwise cause the singularity, this 
fast growth prevents it. Equating these two characteris- 
tic radii gives Tthr ^ j^- As is increased, the region 
where the skewness would cause the singularity moves to 
a smaller radius - into the exponentially-growing part of 
the radial dynamics, and caustics appear. This is the 
picture of the singularity formation in a monostable sys- 
tem: skewness, rather than slowing-down of paths, is the 
key property responsible for the birth of singularities - 
as long as fiuctuational paths move slowly-enough in the 
relevant part of the phase space. Due to the fact that 
our system is minimal - it contains damping and nonlin- 
earity and has only three control parameters, we believe 
this picture may be generic. 

There is a unique correspondence between singulari- 
ties in the pattern of optimal paths and the singularities 
in the probability distribution [17 , as mentioned in the 
Introduction. Indeed, it has been pointed out almost 
twenty years ago by Graham and Tel 31 that if there 
exists a nonequilibrium potential, it will be generically 
non-smooth. The skewness that is the key to our picture 
is a consequence of the loss of integrability - at finite T 
the energy of the auxiliary dynamics is the only integral 
of motion, while at = the angular momentum is 
also conserved. One may ask if any loss of integrability 
will generically lead to singularities? This work demon- 
strates that a loss of integrability is not always sufficient 
for formation of singularities - even when the auxiliary 
dynamics is non-integrable, the nonequilibrium potential 
is everywhere smooth if the driving is not too strong. The 
threshold value depends on the speed of the mean-field 
dynamics at equilibrium. 

The results of this work have been obtained with M. 
I. Dykman, and will be submitted to a peer-reviewed 
journal jointly. The author thanks Michael Cross and 
Michael Khasin for useful discussions; Carl Goodrich, 
Tony Lee and Eyal Kenig for very helpful proof-reading 
suggestions, and the NSF grant number DMR-0314069 
for support in the early stages of this work. 
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Appendix A: Dimensionless variables P —P. The rest of the dimensionless variables and 

parameters are related to the physical ones as follows: 

Variables (Q,P) are related to {q^q) via a canonical 7 
transformation: ^ ~ ^^Jo 7^ (^'^/ 



2\n-ujo\ 

V3 1 



/a 



D = T-TT^ Z2^kT (A5) 



a/ I a/ [P cos m - Q sin m] (A2) Scg (fi - c^o) 



T = -cjol^ = 7^ (A6) 
If a complex convention is used (as in ^2): q = ue^^^^ + 7Cc;q 
c.c. and q = iVtue^^^ -\-c.c., the resulting amplitude equa- = ^ ^^^^ 

tions are related to the present ones by the relationship 
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